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ABSTRACT 


This report contains a brief description of the modifications implemented in the 
PAFAC finite element program for simulation of progressive failure in fibrous composite 
materials and structures. Details of the memory allocation, input data, and the new 
subroutines are given. Also, built-in failure criteria for homogeneous and fibrous 
composite materials are described. 
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1. INTRODUCTION 


It is well known that initiation of local failure in fibrous composites may occur at an 
overall load which is only a fraction of the ultimate strength. For example, Johnson et. al. 
(1983) detected first fiber failure in center notched 0° BfAl composite specimens at about 
60% of the ultimate load. Furthermore, cracks may be caused in the matrix by thermal 
stresses generated during fabrication or after sustained cyclic plastic straining (Dvorak and 
Johnson 1980). Such localized damage states cause reduction of the composites overall 
stiffness, but will not significantly affect the overall strength of a given structure until large 
scale damage develop in the system. The failure process is therefore progressive and must 
be considered in any realistic model. In particular, finite element procedures which are 
capable of performing progressive failure analyses are very useful in predicting failure of 
structures with complex geometry and nonlinear constitutive behavior. 

Several years ago, the finite element program PAC78 (Bahei— El— Din 1979, 

Bahei— El— Din et al. 1981) was developed to perform elastic— plastic analysis of fibrous 
composite structures. Constitutive equations of a unidirectionally reinforced composite 
material were derived from the vanishing fiber diameter model (Bahei— El— Din 1979, 
Dvorak and Bahei-El— Din 1982). The program was later modified under NASA Langley 
sponsorship to perform plastic and failure analysis of composites. This resulted in the 
PAFAC program (Bigelow and Bahei— El— Din 1983) which is currently in use at NASA 
Langley and other government research laboratories and centers. The primary objective at 
that time was to implement a failure criterion based on the fiber stresses in the finite 
element routine and detect the first fiber failure in a given structure under mechanical 
loads. The PAFAC program has been used at NASA Langley for stress analysis of various 
problems in metal matrix composites and laminates (Bigelow 1989, Johnson and Bigelow 
1989, Johnson et al. 1983). 
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The objective of the current project is to develop and implement in the PAFAC 
program the logic required to perform progressive nonlinear/failure analysis for fibrous 
media. The changes in the old PAFAC program, and the failure criteria used in the 
analysis are described in this report. First, we give an overview of the PAFAC program. 
Next, the logic implemented in the program to simulate failure of elements is described. 
Finally, the memory allocation, input data, and the new user— defined subroutines added to 
the program are given. 


2. OVERVIEW OF THE PAFAC FINITE ELEMENT PROGRAM 

PAFAC is a general purpose three-dimensional finite element program. It can be 
used for linear elastic analysis or nonlinear elastic— plastic analysis. Two material types are 
acceptable to the program, homogeneous isotropic materials, and fiber— reinforced 
composite materials. The composite material model is a three-dimensional continuum 
model with unidirectional elastic fibers and an elastic— plastic matrix. The matrix could be 
nonhardening or exhibiting kinematic hardening behavior. The homogeneous material 
model is a three-dimensional continuum model exhibiting kinematic hardening behavior. 
For either material model, the Mises yield criterion and the hardening rule of Prager with 
Ziegler’s modification are employed. A restart capability is built in the program which 
allows analysis of the problems with previous loading history. Loads handles by the 
PAFAC program are proportional mechanical loads with any loading, unloading, and 
reloading sequence. 

The PAFAC program does not include any pre— or post— processors. The users, 
however, may prepare their own subroutines for automatic generation of nodal coordinates, 
element properties and connectivity, boundary conditions, and yield and hardening 
information for all material types. 
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The three-dimensional solid element (eight-node hexahedron) is the only type that 
the PAFAC program uses to span the material volume. By means of the deflection 
boundary conditions, any type of deformation may be simulated, for example, beams, 
plates, plane strain, etc. Governing equations for the finite element solution are generated 
in the context of the displacement method of analysis. Linear equilibrium equations are 
solved by Cholesky method. Nonlinearities are handled by a Newton— Raphson type 
iterative scheme. In its present form, the PAFAC program predicts initial failure loads in 
fibrous composites using a quadratic failure criterion for the fiber phase function of the 
axial stress and the longitudinal shear stress. 

Description of the input data and the output of the PAFAC program can be found 
in the user’s manual prepared by Bigelow and Bahei— El— Din [1983]. Solutions for a 
variety of laminated plate problems can be found in Bigelow (1989), Johnson and Bigelow 
(1989), and Johnson et al. (1983). 


3. DESCRIPTION OF THE MODIFIED PAFAC PROGRAM 
3.1 Progressive Failure Analysis 

In the modified PAFAC program, element failure is examined upon completion of a 
loading step. The user can choose to examine failure using failure criteria defined by the 
program, or specify the failure criteria by supplying two subroutines, UFAIL to describe 
the failure criteria and check element failure, and DFAIL to provide the input data 
required by the user— defined failure criterion. The form of the failure criteria specified by 
the program depends on the type of material. For isotropic materials, the following failure 
envelope is specified in the stress space: 
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the fiber stress. In either case, only the axial normal stress and the longitudinal shear 
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Here, a i j is the overall stress and cry is the fiber stress specified in a Cartesian coordinate 
system x k , k = 1,2,3, such that the x 3 — axis coincides with the fiber axial direction, and 
coincides with the transverse plane. The symbols cr u , a u , denote overall strength of the 
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fibrous composite under axial normal stress and longitudinal shear stress, respectively. 
Similarly, a uf , r uf denote the fiber strength under axial normal stress and longitudinal 
shear stress, respectively. 

When element failure is detected by the program, the forces equilibrating the 
stresses supported by the failed element are retrieved from the common array, and their 
norm, RNORM, is computed as the sum of the absolute value of the individual 
components. This force vector is re— applied to a new mesh in which the failed element is 
removed, and a new overall stiffness matrix is formed. The force vector, however, is 
applied in small increments, such that the norm of each increment does not exceed the 
norm, PNORMB, of the overall load increment applied to the structure prior to the 
element failure. Specifically, the program scales down the residual load vector computed at 
element failure by the ratio PNORMB/RNORM if it is less than unity, otherwise, the 
residual load vector is applied back to the modified mesh as one load increment. During 
redistribution of the residual load vector found for a failed element, the overall loads are 
held constant. If in this process more elements fail, the residual load vectors corresponding 
to the stresses in the failed elements are computed and the overall residual load vector is 
updated. After redistribution of the total stress supported by the failed elements, a new 
overall load increment is applied and the process for evaluation of the failure criterion in 
the element and redistribution of the stresses, if necessary, is repeated. 

3.2 Modifications to Memory Allocation and Input Data 

Minor modifications are made in the memory allocation and input data. Table 1 
shows the common map of the general constants and arrays in the common block of the 
modified PAFAC program. The input data for a standard job are given in Tables 2-4. 

For a restart job, the input data are given in Tables 5,6. 
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3.3 New User— Defined Subroutines 

Two user— defined subroutines are added to the new PAFAC program, UFAIL and 
DFAIL, which specify the failure criterion and read the failure parameters. Listing of the 
DFAIL subroutine is given in Fig. 1, and listing of the UFAIL subroutine is given in Fig. 2 
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Table 1 Common map of the general constants and arrays in the common block 

of the modified PAFAC program 


Location 

in 

Common 

Symbol 

Description 

1 

IN 

Total number of nodal points 

2 

IBN 

Total number of input units for deflection boundary 
conditions (dbc) 

3 

IT 

Total number of elements 

4 

IP 

Total number of concentrated loads 

5 

IPRS 

Number of pressure units 

6 

ITYPE 

Indicator for material type 

7 

IMAT 

Number of material types 

8 

INX 

Number of last link to be executed 

9 

IMFI 

Number of fiber angle types 

10-17 

Ni 

Labels of current element vertices 

18 

M 

Label of current element 

19 

ITY 

Common size in words 

20 

ITEM 

Temperature change type number 

Material type number 

Number of reduced equilibrium equations 

21 

IMET 

22 

ISUM 

23 

ISHUF 

Node relabelling indicator 

24 

10 RD 

Number of words for reduced stiffness matrix 

25 

IORD1 

IORD+1 

26 

ACEL 

Body force per unit volume 

27 

INP 

Indicator for output level 

28 

IPBG 

Integer constant for element load vector 

29 

IPEN 

Integer constant for element load vector 

30 

CONS 

Constant for element load vector 

31 

IU 

Pointer for diagonal element count array 

32-34 

G1,G2,G3 

Direction cosines of acceleration vector 

35—42 

43 

J1,J2,..J8 

Pointer for J1W, J2W,N1,N2,..N6 arrays 
not used 

44 

IBB 

Pointer for IBB array 

45 

IBO 

Pointer for IBO array 

46 

IID 

Pointer for material constant array 

47 

IIA 

Pointer for thermal coefficient array 

48 

IDT 

Pointer for temperature change array 

49 

ICFI 

Pointer for fiber angle array 

50 

IXX 

Pointer for x— coordinate array 

51 

IYY 

Pointer for y— coordinate array 

52 

IZZ 

Pointer for z-coordinate array 

53 

IIC 

Pointer for constant array of dbc units 

54 

IDEF 

Pointer for unknown deflection array 

55 

1ST 

Pointer for overall reduced stiffness matrix 

56 

IIS 

Pointer for element stiffness matrix 

57 

IERR 

Error indicator 

58 

DT 

Value of temperature change of element 

59-£l 

AL1.AL2, 

AL3 

Coefficients of thermal expansion of element referred to 
material axes 
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Table 1 Common map of the general constants and arrays in the common block 
of the modified PAFAC program (Cont.) 


Location 

in 

Common 

Symbol 

Description 

62-82 

D21i 

Material matrix of an element 

83-106 

Pi 

Load vector of an element 

107-130 

UVi 

Deflections of an element due to temperature change 

131-138 

Xi 

Overall x-coordinates of element’s vertices 

139-146 

Yi 

Overall y— coordinates of element’s vertices 

147-154 

Zi 

Overall z-coordinates of element’s vertices 

155-161 

XDi 

Xi-X t 

162-168 

YDi 

Yi-Yi 

169-175 

ZDi 

Zi-Z, 

176 

IMES 

Indicator for element descriptors input 

177 

IBUN 

Indicator for boundary conditions input 

178 

ICOR 

Indicator for nodal coordinates input 

179 

IPIR 

Unit number for restart jobs input data 

180 

PRES 

Pressure value for an element 

181 

IPR 

Pointer for pressure array 

182 

ITAS 

Unit number for intermediate storage 

183 

JMFI 

Fiber angle type number 

184 

JMMZ 

Pointer for dbc constant array 

185 

JMMY 

Pointer for dbc second pair array 

186 

JMMX 

Pointer for dbc first pair array 

187 

JPRS 

Pressure type number 

188 

J10 

Pointer for Ns array 

189 

J9 

Pointer for N 7 array 

190 

191-199 

200-325 

ISDT 

Number of temperature change types 
General usage area 
User’s area 

326 

IO 

Formatted input unit number 

327 

MO 

Formatted output unit number 

328 

MAXI 

Substep printout indicator 

329 

MAXSI 

Maximum number of substeps 

330 

NPRIN 

Printout indicator for stress link 

331 

NSTAV 

Strain computation method indicator 

332 

ER1 

Convergence criterion for substep iteration 

333 

KHARD 

Strain— hardening indicator 

334 

IHARD 

Indicator for hardening information input 

335 

NAA 

Count of vectors for saving element’s response 

336 

LHD 

Pointer for material plastic properties array 

337 

IELAS 

Indicator: Oselastic element; l=otherwise 

338-358 

LLLi 

Pointer of NAA arrays 

359 

ITER 

Step counter 

360 

ITERS 

Substep counter 

361 

PNORMB 

Norm of overall load vector 

362 

PNORM 

Norm of residual load vector 

363 

IDIV 

Indicator for subdivision of hexahedrons 

365 

CUMS 

Cumulative load factor of current step 

366 

DESC 

Incremental load factor of current step 


12 






Table 1 Common map of the general constants and arrays in the common block 
of the modified PAFAC program (Cont.) 


Location 

in 

Common 

Symbol 

Description 

367 

LLD 

Pointer for load factor-step number array 

368 

LCDP 

Number of data points in LLD array 

369 

MPR 

Internal indicator for stress link printout 


JJ 

Indicator: 0=elastic— plastic; Inelastic 

371 

IPOT 

Unit number of plotting file 

372 

YJ2 

Effective stress 


H 

Plastic tangent modulus 

374 

ITAX 

Unit number for cumulative data 

375 

ZZ 

Scale factor of the step 

376 

IAMB 

Ambient temperature of the structure 

377 

IYIT 

Indicator for isothermal materials 

378 

PSCAL 

Storage for original scale factor 

379 

ISECT 

Internal indicator for isothermal processes 

380 

IYEL 

Indicator for generation of stiffness matrix 

381 

KTERS 

Internal storage for ITERS 

382 

QNORMB 

Internal storage for PNORMB 

383 

D 

Internal storage 

384 

IDD 

Job identification 

385 

IRR 

Job counter 

386 

NI 

Count of elements in plastic zone 

387 

ICOUNT 

Internal storage 

388 

ISCAL 

Indicator for scaling in a restart job 

389 

IFLAG* 

Indicator: 0=regular load increment; >0 failure increment 


IFAL* 

Failure criterion indicator 

391 

NOFAL* 

Pointer for element failure indicators 

392 

IIF* 

Pointer for parameters of failure criterion 


* New item. 
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Table 2 Input items of a standard PAC91 job 


Input 

Item 

No. 

Conditions 

Input Item Contents 

Format 

1 


(Bi, i=l,20), This line contains 
any alphanumeric message except 
the first four columns that 

should not contain the literal 
.????! 

20A4 

2 


IN ,IT ,ITYP E ,IBN ,IP ,IP RS,IM AT, 
ISDT,IMFI,LCDP,INX,INP,ISHUF, 
ICOR,IBUN,IMES,IPIR,KHARD,IHARD, 
mT,ITAS,IAMB,NAA,IFALGl,G2,G3, 
ACEL (See Table 3 for details) 

214,11,314,412, 
1111, 15, 212, 6x, 
3F5.4.E10.3 

3 

itype=o 

(i,Ei,i/i,ai, i=l,IMAT) 

(3(12, 3E8. 2)) 


ITYPE=1 

i=l,IMAT) 

(I2,F5.3,4E12.5) 

4 

KHARD=1 

1<IHARD<9 

(i,j,(fk,^k, k=l,j), i=l,IMAT), 
2<J<10*IHARD 

(2I3/(10E8.3)) 


KHARD=1 

IHARD=0 

Input should be prepared as 
required by user’s subroutines 
YHARP,YHAR2P,YHAR4P 



KHARD=0 

No input 


5 


(Ki,Ci, i=l,LCDP) 

(8(12, E8.2)) 

6 

1<IPRS<99 

(i,Pi, i=l,IPRS) 

(8(12, E8.2)) 


IPRS=0 

No input 


7 

1<IMFI<99 

i=l,IMFI) (8(I2,E8.5)) 

(see Bigelow and Bahei-El— Din (1983) for details) 


IMFI=0 

No input 


8 

ICOR=0 

(j>Xj,yj,Zj, j=l,IN) 

(2(I4,3E12.4)) 


ICOR=l 

Input should be prepared as 
required by user’s subroutine 
CORG 
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Table 2 Input items of a standard PAC91 job (cont.) 


Input 

Item 

No. 

Conditions 

Input Item Contents 

Format 

9 

IBUN=0 

(ik)jk}ik,jk,&k, k=l,IBN) 

(5(14,11,14 

I1.F6.0)) 


IBUN=1 

Input should be prepared as 
required by user’s subroutine 
BUNG 


10 

IMES=0 

(MM in ,JlW in) J2Wni,Nl....N8, m=l,IT) 
MMin = — (m—(m/ 1000)* 1000) interpreted 
in Fortran integer arithmetic 
SGHS6’ 

J1W=100*IMET+JMFI, 

J2W=100*JPRS+ITEM 

(2014) 


IMES=1 

Input should be prepared as 
required by user’s subroutine 
MESG 


11 

ISHUF=0,1 

No input 



ISHUF=2 

(Ni, i=l,N) 

(2014) 

- 

ISHUF=3 

(Ni,IMAXi, i=l,IN) 

(2014) 

12 

1<IP<9999 

(il,ji,Pi, 1=1, IP) 

(5(14,11 ,E1 1.4)) 


IP=0 

No input 


13* 

IFAL=0 

Input should be prepared as 
required by user’s subroutine 
DFAIL 



IFAL=— 1 

(n,(ri, i=l,IMAT) 

(2E10.3) 


IFAL=1 

ITYPE=1 

(n f ,<7i f , i=l,IMAT) 

(2E10.3) 

14 


MAXI,MAXSI,NPRIN,NSTAV,ER1,IDIV, 
IPOT,ITAX,IYEL,IDD,IRR,TEST 
(See Table 4 for details) 

412, E12.5, 11,12 

I1,I4,A4,I2,43X, 

’END’ 


* New item. 


15 










Table 3 Description of input item 2 of a standard PAC91 job 


Name Columns Format Range Description 


IN 

1-4 

14 

8-9999 

IT 

5-8 

14 

1-9999 

ITYPE 

9 

11 

0-1 


IBN 

10-13 

14 

3-9999 

IP 

14-17 

14 

0-9999 

IPRS 

18-21 

14 

0-9999 

IMAT 

22-23 

12 

1-99 

ISDT 

24-25 

12 

0 

IMFI 

26-27 

12 

0-99 


LCDP 

28-29 

12 

1-25 

INX 

30 

11 

1—4 

INP 

31 

11 

0-2 

ISHUF 

32 

11 

0-3 


ICOR 

33 

11 

0-1 

IB UN 

34 

11 

0-1 

IMES 

35 

11 

0-1 

IPIR 

36 

11 

0-9 

KHARD 

37 

11 

0-1 

IHARD 

38 

11 

0-9 


Total number of mesh points 
Total number of elements 
Material indicator: 

0 — homogeneous, isotropic, 

1 - fibrous composite 

Total number of deflection boundary 
condition units 

Total number of concentrated load units 
Total number of different pressures 
Total number of different materials 
Total number of different temperature 
changes 

Total number of different fiber angle 
types: 

0 if ITYPE=0; 

1-99 if ITYPE=1 

Total number of load factor-step number 
points 

Number of last link to be executed 
Printout indicator: 

0 — minimum output; 

1 — intermediate output; 

2 — detailed output 
Relabelling indicator: 

0 — no relabelling; 

1 — iterate to relabel without reading 

data; 

2 — read data and iterate to relabel; 

3 — relabel as shown on data lines 
Indicator for nodal coordinate generation: 

0 — read from data lines; 

1 - subroutine CORG 
Indicator for deflection boundary 
conditions generation: 

0 — read from data lines; 

1 — subroutine BUNG 
Indicator for element descriptors 
generation: 

0 — read from data lines; 

1 — subroutine MESG 

Unit number for storage of information 
necessary for restart jobs 
Strain hardening indicator: 

0 — nonhardening (only if ITYPE = 1); 

1 — kinematic hardening 
Indicator for hardening information: 

1—9, piecewise linear; 

0 , not piecewise linear or nonhardening 
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Table 3 Description of input item 2 of a standard PAC91 job (cont.) 


Name 

Columns 

Format 

Range 

Description 

IYIT 

39 

11 

0 

Indicator for yield stress— temperature 
information 

ITAS 

40 

11 

1-9 

Unit number for intermediate storage 

IAMB 

41—45 

15 

0 

Ambient temperature of material 

NAA 

46—47 

12 

15,21 

Number of arrays for saving element 
information: 





15 if ITYPE=0; 
21 if ITYPE=1 

IFAL* 

48—49 

12 

-1,0,1 

Failure criterion indicator: 

—1 : program-defined based on overall 





stress 





0 : user— defined via subroutine UFAIL 





1 : program-defined based on fiber stress 
(only if ITYPE=1) 

G1 

56-60 

F5.4 

-1 to 1 

X— component of unit acceleration vector 

G2 

61-65 

F5.4 

-1 to 1 

Y-component of unit acceleration vector 

G3 

66-70 

F5.4 

-1 to 1 

Z— component of unit acceleration vector 

ACEL 

71-80 

E10.3 

any 

Magnitude of acceleration vector times 
unit mass, i.e. unit weight 


* New item. 
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Table 4 Description of input item 14 of a standard PAC91 job 
and input item 6 of a restart PAC91 job 


Name Columns Format Range Description 


MAXI 

1-2 

12 

0-99 

MAXSI 

3-4 

12 

-9 to 99 

NPRIN 

5-fi 

12 

0-2 


NSTAV 

7-8 

12 

0-1 

ER1 

9-20 

E12.5 

any > 0 

IDIV 

21 

11 

0-1 

IPOT 

22-23 

12 

0-9 

ITAX 

24 

11 

1-9 

IYEL 

25-28 

14 

1-9999 

IDD 

29-32 

A4 

any 

IRR 

33-34 

12 

1-99 

TEST 

78-80 

’END’ 



Printout indicator for stress link: 

> 0 — printout is produced at every step 

and at every substep; 

0 — printout is produced at the steps 
given in input item 5 and at the 
end of failure steps 

Maximum number of substeps; if < 0, 10 
substeps are assumed 
Printout level for stress link: 

0 — stresses and forces are printed 

1 — stresses, forces, and strains are 

printed 

2 — stresses, forces, strains and 

displacements are printed 
Indicator for strain computation from 
displacements: 

0 — best fit strain tensor is obtained 

1 — computations by user’s subroutine 

STNA 

Convergence tolerance; if 0, 0.01 is 
assumed 

Indicator for subdividing hexahedrons: 

0 — in two ways; 

1 - in one way 

Unit number for storage of information 

necessary for plotting 

Unit number for storage of cumulative 

results 

Number of elements allowed to yield in a 
substep without generation of stiffness 
matrix 

Identification; not blank if IPOT > 0 
Sequential job label: 

1 for standard job 

> 1 for restart job 
Indicator for end of current job 
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Table 5 Input items of a restart PAC91 job 


Input 

Item 

No. 

Conditions 

Input Item Contents 

Format 

1 


(Bi, i=l,20), This line contains 
any alphanumeric message except 
the first four columns that 
should contain the literal *????* 

20A4 

2 


NIP,NIPRS,NINP,NIHARD,JTAS,NIPIR 

ISCAL,NLCDP 

(See Table 6 for details) 

14, 12, 10X, 511 ,12 

3 


(Ki.Ci, i=l,LCDP) 

(8(12, £8.2)) 

4 

NIPRS=IPRS*, 

1<IPRS<9999 

(i, Pi , i=l,NIPRS) 

(8(12, E8.2)) 


NIPRS^IPRS 
or IPRS=0 

No input 


5 

NIP=IP* 

1<IP<9999 

(il.jl,Pl, 1=1, NIP) 

(5(I4,I1,E11.4)) 

- 

NIP^IP 
or IP=0 

No input 


6 


MAXI,MAXSI,NPRIN,NSTAV,ER1,IDIV, 
IPOT ) ITAX,IYEL,IDD,IRR,TEST 
(See Table 4 for details) 

4I2,E12.5,I1,I2 

I1,I4,A4,I2,43X, 

’END’ 


* IPRS = Number of pressure types of the preceding job 

* IP = Number of concentrated load units of the preceding job 
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Table 6 Description of input item 2 of a restart PAC91 job 


Name 

Columns 

Format 

Range 

Description 

NIP 

1-1 

14 

0-9999 

Total number of concentrated load units 

NIPRS 

5-6 

12 

0-99 

Total number of different pressures 

NINP 

17 

11 

0-2 

Printout indicator: 

0 — minimum output; 

1 — intermediate output; 

2 — detailed output 

NIHARD 

18 

11 

0-9 

Indicator for hardening information: 

1—9, Information of standard job is used; 
0 , user’s subroutines YHARP, YHAR2P, 
YHAR4P are used 

JTAS 

19 

11 

1-9 

Unit number for restart input information 

NIPIR 

20 

11 

0-9 

Unit number for storage of information 
necessary for succeeding restart job 

ISCAL 

21 

11 

0-1 

Scaling indicator: 

0 — load factors of input item 3 are not 

scaled with SCAL* before their use; 

1 — load factors of input item 3 are scaled 

with SCAL before their use 

NLCDP 

22-23 

12 

1-LCDPt 

Total number of load factor-step number 
points 


* SCAL = Load factor for initial yielding computed in the initial standard job 

* LCDP = Number of data points of load factor-step number in the initial standard job 
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c 


SUBROUTINE DFAIL 

To read in user— provided input data related to failure 


IMPLICIT DOUBLE PRECISION (A-H,0-Z) 
COMMON/ONE/IN 
COMMON/TWO/XIN 
DIMENSION IA(401),AA(401),S(1),IS(1) 


EQUIVALENCE 

EQUIVALENCE 

EQUIVALENCE 


(IN,IA 

IA(7) 




IO),(AA(392),IIF) 


IA(401),S(l)),(AA(40l),IS(l)) 


C..The user provides here statements which read failure parameters for every material 
C..type up to IMAT. The data are to be stored in array AA in the area specified by the 
C.. pointer IIF such that all parameters for material type 1 are listed in the first 10 
C.. positions followed by all parameters of material type 2 in the next 10 positions, etc. 
C..The read statement should specify the input unit number as 10. 


C.. Example, read one material parameter 

DO 10 1=1, IMAT 
IPOINT=IIF+(IMAT-1)*10 
READ (10,*) AA(IPOINT+l) 

10 CONTINUE 

C..End of Example 


RETURN 

END 


Fig. 1 Listing of user— defined subroutine DFAIL 
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SUBROUTINE UFAIL(KFAIL) 


To provide the user’s failure criterion and check if it is satisfied for element M 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

COMMON/ONE/IN 

COMMON/TWO/XIN 

DIMENSION IA(401),AA(401),IS(1),S(1) 

DIMENSION LLL(21) 

EQUIVALENCE (IN,IA(1)),(XIN.AA(1)) 

EQUIVALENCE (lA( 6 ),ITYPE),(lA( 18 ),M)^IA( 21 ),IMET),(IA( 46 ),nD) 


EQUIVALENCE (IA(338),LLL(l)) > (AA(392) > nF) > (IA(401),IS(l)),(AA(401) > S(l)) 
EQUIVALENCE (AA(236),DIR(l,l)) l (AA(274),E66(l,l)) 

DIMENSION DIR(3.3),E66(6,6) 

DIMENSION SSO(6),SSL(6),SSM(6),SSF(6) 


KFAIL=— 1 


C.. Overall stresses of the current element are stored in array SSO. For a composite 
C.. element, the stresses contained in SSO are described in the global coordinate system 
C. .while those given in array SSL are described in the local coordinate system where 
C. .direction 3 is the fiber axial direction, and plane 1—2 is the transverse plane. The 
C. .corresponding matrix and fiber stresses are stored in arrays SSM and SSF, 

C. .respectively. The stress components are listed in arrays SSO, SSL, SSM, SSF in the 
C.. order S11,S22,S33,S12,S13,S23. 


DO 10 1=1,6 
SSO(I)=AA(LLL(I)+M) 
10 CONTINUE 


30 


2 


IF (ITYPE.GT.O) THEN 
CALL CODI4 
CALL TRANS4 
DO 30 1=1,6 
SSL(I)=0.0D0 
DO 30 J=l,6 

SSL(I)=SSL(I)+E66(I,J)*SSO(J) 
CONTINUE 
IIDI=IID+(IMET— 1)*5 
VF=AA(IIDI+1) 

VM=1.0D0— VF 
DO 20 1=1,6 

SSM(I)=AA(LLL(I+9)+M) 

SSF(I)=(SSL(I)-VM*SSM(I))/VF 

CONTINUE 

ENDIF 


Fig. 2 Listing of user— defined subroutine UFAIL 
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C..This part is prepared by the user to check the failure criterion for the current 
C..element. Input data entered by the user in subroutine DFAIL are stored in array AA 
C..in the area specified by the pointer IIF in the following order: parameters for 
C.. material 1 are located in the first 10 positions, and those for material 2 are located in 
C..the next 10 positions, etc. If the failure criterion is satisfied, set KFAIL=1. 


C.. Example, failure occurs when the axial stress in the fiber reaches a uf 

IPOINT=IIF+(IMET— 1)*10 
SULT=AA(IPOINT+l) 

IF(DABS(SSF(3)).GE.SULT) KFAIL=1 

C..End of Example 

RETURN 

END 


Fig. 2 Listing of user— defined subroutine UFAIL (Cont.) 
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